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^ ; Abstract 

CN ' We report and explain a convective phenomenon observed in molecular dy- 

K.^ namics simulations that cannot be classified either as a hydrodynamics in- 

QQ ' stability nor as a macroscopically forced convection. Two complementary 

o' 

arguments show that the velocity field by a thermalizing wall is proportional 
to the ratio between the heat fiux and the pressure. This prediction is quan- 
C^ ' titatively corroborated by our simulations. 
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Free thermal convection — driven by buoyancy or by surface tension — is a perfectly well 
understood phenomenon derivable from Navier Stokes equations [|1],0]. Simulations of free 
thermal convection by means of Molecular Dynamics (MD) techniques can be achieved with 
systems with as few as 10^ particles and already these small systems exhibit hydrodynamic 
behavior as seen for example in |^- ^. Moreover MD is useful to study fluid phenomena 
at the microscopic level without having to make assumptions concealed behind the Navier- 
Stokes equations such as the Fourier law, Newton's law of viscosity and local thermodynamic 
equilibrium. 

In the following we are going to present a unusual convective phenomenon (not pre- 
dictable by Navier-Stokes equations) related to the variation of the temperature field in 
one mean free path i through the adimensional parameter iVT/T. This parameter can 
be interpreted as a measure of how far from local thermal equilibrium the system is at a 
given point. When effects violating the local thermodynamic equilibrium are present one 
should question the very concept of temperature but we will manage without engaging in 
such delicate matters. 

The convective phenomenon that we are reporting takes place when there is a temper- 
ature gradient parallel to a thermal wall. The mechanism can be sketched as follows: the 
particles that approach a point P of the wall come from an anisotropic distribution while 
the particles that hit the wall at P come back to the system with a distribution which is 
isotropic, or at least less anisotropic than the incoming flux. A careful assessment of the 
difference between the incoming and outgoing fluxes at P yields the conclusion that there is 
a net mass flux parallel to the wall. We have observed this phenomenon in MD simulations 
and have made a theoretical estimation of its value. In real experiments the effect will be 
small but it should be observable in a rarefied gas. 

We have made MD simulations of a two dimensional gas of hard disks in a square 
box, using our own efficient algorithm ^ and the carefully devised measurement routines 
described in 0. Each numerical experience consisted of two runs: (1) The system with 
periodic vertical walls was subjected to a temperature difference, relaxed for 200 thermal 



diffusion times tr, and then the temperature profile T{y) was carefully measured for another 
200 tj-- (2) A second simulation was run under the same conditions as in (1), except 
that a periodic permeable thermal vertical wall was added. This new wall was defined to 
have the previously obtained profile T{y), namely each particle hitting the vertical wall 
emerged on the other side of the box with a velocity taken from a heat bath at the local 
temperature T{y). The sign of the vertical component of this velocity was random, and 
therefore microscopically it is a non slip boundary condition, in the sense that the emerging 
particles do not remember the velocity with which they came. Letting the particles pass 
through the (periodic) wall is totally irrelevant to the resulting phenomenon, but it helps 
reducing the boundary effects near the wall. Again the system was relaxed for 200 t^, and 
then measurements were averaged in time during the next 600 ty. The measurements where 
done dividing the system in square cells. Densities and the velocity field were measured in 
every cell, and fluxes were measured across the cell walls. 

Units are chosen such that particle's mass and diameter, the Boltzmann constant, and 
the temperature at the bottom, m, D, ks and T^ respectively, are fixed to unity. With 
this particular choice of units the lengths are in diameter units, the temperature in energy 



units, and the time in units of JmD'^/kBTi,. The control parameters of each simulation are 
the number of particles N, the bulk number density hb, and the temperature at the top Tt 
where Tt<Ti,. 

Our main simulation considered a system of A^ = 1444 hard disks, bulk number density 
ub = 0.05, implying a box side of 170 and a mean free path of about 7, and at the top the 
temperature was fixed to be Tt = 0.1. 

The main observation is the following: a convective current stabilizes in the neighborhood 
of the vertical wall, moving towards the warmer zone. In figures |1] and |^ it is possible to see 
the velocity field v and the mass flux mnv. At the bottom the convective current necessarily 
bends towards the center to come up along the central part of the box. Since the gas is 
highly compressible, the eye of the convective rolls are far from the expanded hotter zone. 
The velocity component Vy in the cells by the vertical wall is almost constant, and its average 



was 



Vy = -0.015 ± 0.002 [observed] (1) 

after excluding 10 cells in the upper and lower extremes, with 76 x 76 the total number of 
cells. 

The same convective phenomenon was observed in all the other situations we simulated: 
(i) n = 0.01, N = 8100, Tb = 1.0, and Tt = 0.01, (ii) n = 0.25, A^ = 1444, Tb = 1.0, 
and Tt = 0.1. The velocity component Vy measured near the vertical walls in (i) was 
Vy = 0.014 ±0.003 and it shows the same behavior as the preceding simulation. The vertical 
component of the velocity in (ii) however is not longer constant, it increases with height. 
The theoretical derivations that we make below are not applicable to this denser case but it 
is interesting to observe that the phenomenon still exists. 

Finally we made another simulation in which the temperature profile of the thermalizing 
wall was not the one obtained from a first run but rather T{y) was chosen arbitrarily to be 
a smooth monotonic profile. In this case we used n = 0.05, A^ = 1444 and Tf = 0.1. Again a 
convective current was created with similar characteristics to the previous ones. This result 
and the theoretical calculations below suggest that to obtain this convective motion, it is 
enough to have a temperature gradient parallel to a thermalizing wall so that each particle 
emerging from the wall comes from a (at least partially) thermalized distribution. 

This kind of convective motion can not be obtained within the usual frame of Navier- 
Stokes hydrodynamics, unless these equations are solved imposing by hand a nonvanishing 
tangential velocity as boundary condition at the vertical wall. That however would be arti- 
ficial because hydrodynamics is built using only the first three momenta of the distribution 
function which are not enough information to take into consideration the phenomena that 
we are reporting. 

In what follows, we give two heuristic and complementary derivations for a rarefied two 
dimensional hard disk gas, one based on local nonequilibrium distribution functions, and 
the other one based on the mean free path theory of transport. Both derivations yield 



essentially the same prediction for the velocity field near the vertical thermalizing wall. 
What the calculations below imply is that the vertical wall exerts an effective tangential 
force on the gas such that a velocity field — proportional to the ratio between the heat fiux 
and the pressure — pointing against the heat fiux is established. 

The basic idea behind the following two derivations is that particles hitting the thermal- 
izing wall at a point P (see figure ^ come from an anisotropic nonequilibrium environment, 
while the particles emerging from P come from an equilibrium isotropic distribution. It is 
understandable then that some fiuxes do not necessarily cancel and in particular we prove 
that there is a net velocity field. 

Nonequilibrium Interpretation: The velocity distribution function near a point y 
of the thermalizing wall (figure |^) has two contributions: (a) one from the particles that 
come towards y from a nonequilibrium velocity distribution and (b) the other one from the 
outgoing particles that come from the thermal bath at y. The nonequilibrium distribution 



function for a system under a heat fiux adapted from |T^ to the case of a two dimensional 
system is 



I m 

/neq = I ^ + ^ 



mv'^ 
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^•g /eq (2) 
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where /eq is the usual Maxwellian distribution. Then the velocity distribution near the right 
wall is 

^ I /neq(^?) v, > 
[ feq{v) v^<0 
Using this distribution, the x-y component of the stress exerted by the fiuid against the 
wall {axy = mn{vxVy) f) is 



Then the force per unit length exerted by the wall on the fiuid is, by the action-reaction 
principle, the negative of the previous expression, hence it points antiparallel to the heat 
flux. 



The velocity near the wall can be estimated using Newton's law {axy = rjdvy/dx where r) 
is the shear viscosity) and the assumption that this velocity decays at distances comparable 
with the mean free path i and it is, 

fa 

Vy = --^ (5) 

T) 

Replacing the expressions for the shear viscosity, the mean free path, and using the 



equation of state for an ideal gas yields |11 



Vy = -g- (6) 

Kinetic Interpretation: Let us consider the mass flux balance at an arbitrary point P 
of the wall, as shown in figure |^. The mass flux coming from an angle between and ip + dcj) 
with respect to the normal to the wall and reaching P is 

It 

dj'=mn\ — a((j)) (cos (f), — sin (f)) d(f) (7) 

where n and T are the number density and temperature at the points where the particles 
come from whereas a(0) is a geometrical factor depending on the incident angle. We do not 
give a value for a{(f)) since it is well known that the mean free path theory of transport is 



too simple to produce the correct numerical factor ||rT|. Since the number density is small 
then p = nT. 

The combined inward mass flux from the directions (p and —0 to P is then 

dv,-n((l)) = \/mpa((j)) cos(7!) , H , x — sind) \ , , y 

The combined outward mass flux from directions and —(f) can only be in the x direction, 
since the flux comes from the local equilibrium distribution at the wall, and it is rfj^^t i^) ~ 



— mn J ^:^ cos (j)xd(j). The density n in rfj^^t ^^^ ^°^ t)e replaced by p/T, because the 
equation of state is not valid at wall points. This number density is unknown and can be 
determined imposing null net mass flux in the x direction. 
The total flux in the y direction is 



'^-y = Jq (^Jin + dfout) ■ y = ^v^ J^ -^ (8) 



with 

CTt/2 



A = a((f)) sin^ (pdcf) 

Jo 



and where we have used that T(0) = T + i sin (p dT / dy with i the mean free path. It must 
be remarked that the previous expression can be written as 

r ^ fidT\ 

where f^.]^ is the thermal velocity. The previous result implies that a mass flux parallel to 
the temperature gradient is induced near the wall, and it is proportional to the adimensional 
parameter |;^, which is a measure of how far from local thermal equilibrium the system is 
at a given point. 

The velocity field near the wall can be estimated dividing this flux by the mass density, 
using the Fourier law, the expressions for the mean free path, the thermal conductivity, and 
the equation of state for an ideal gas [|r^. The outcome is 

A [Wq , , 

This result predicts the same behavior as (|^). 

From these two heuristic arguments we have shown that the mechanism for this type of 
convection is already present in nonequilibrium thermodynamics and in fact both formula- 
tions, mean free path theory and nonequilibrium local distribution functions, give essentially 
the same result: the velocity field near the wall is proportional to the heat flux, and it is 
independent of the point P — due to the absence of external forces and energy sources, p 
and q are uniform — as it can be appreciated in figure |I]. Our predicted value for Vy from 
our observations of q and p is 

Vy = -0.016 ± 0.003 [predicted] (11) 



which should be compared with (|l|), and Vy = —0.020 ± 0.001 for the simulation with 
n = 0.01, N = 8100, Tfe = 1.0, and Tt = 0.01. 

The extension to three dimensions is straightforward giving essentially the same result, 
indicating that this convective motion could be observed in a rarefied gas. For example using 
gaseous Helium at atmospheric pressure with a temperature gradient of VT = 100[-ft'/c?7i], 
the velocity near the wall that we are predicting is Vy = 3.8[mm/s], but it may be somewhat 
smaller because, in a real experiment, the average flux coming out from every point at the 
thermalizing wall is not totally isotropic. 

Regarding the heat flux q, there has been an interesting recent proposal applicable to 
the case of systems under a heat flux, saying that one should observe a departure from the 
standard Fourier law [|1^] (which actually motivated our series of simulations.) What we 
observe (see figure H) is that the energy flux is consistent with the Fourier law plus the 
kinetic flux \{mnv'^ v) of an ideal gas. The effect predicted in [Q for our system is about 



an order of magnitude smaller that the total flux, and since fluxes are noisier than densities, 
we can not yet see if such an effect exists. 

In summary, we have observed and justified that there exists a free thermal convection, 
which stems from non local effects due to the presence of nonequilibrium distributions. The 
velocity by the wall is proportional to I'WT /T, which is assumed to vanish in hydrodynamics 
(local thermodynamic equilibrium). 
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FIGURES 
FIG. 1. Velocity field measured using MD simulations. The horizontal walls are kept at uniform 

temperature, the warmer wall is at the bottom, and the vertical walls have a different temperature 

at each point. The number of particles, the bulk number density, and the top temperature were 

N = 1444, n = 0.05, and T( = 0.1 respectively. 

FIG. 2. Mass flux measured in the same simulation shown in the previous figure. 

FIG. 3. Contributions to the mass flux from different directions. Particles come from regions 
at different temperatures and emerge with velocities from an isotropic distribution at temperature 
T. 

FIG. 4. Energy flux measured in the same simulation as in previous figures. The x component 
has been amplified four times to show how it is distorted by the convective current. 
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